
#ifndef LBM2D_H
#define LBM2D_H
 
#include<cmath>
#include<cstdlib>
#include<iostream>
#include<vector>

#define Q_DIM 9

using namespace std;

template <typename T> struct vec {
	T x;
	T y;

	double norm() {
		return x*x + y*y;
	}

	T operator*(const vec<T>& vec2) const {
		return x*vec2.x + y*vec2.y;
	}
};

const vec<double> eVector[Q_DIM] = {{0, 0}, {1.,0}, {0,1.}, {-1.,0}, {0,-1.}, {1.,1.}, {-1.,1.}, {-1.,-1.}, {1.,-1.}};
const int invertedIndices[Q_DIM] = {0, 3, 4, 1, 2, 7, 8, 5, 6};
const double wCoeff[Q_DIM] = {4./9., 1./9., 1./9., 1./9., 1./9., 1./36., 1./36., 1./36., 1./36.};

enum direction {
	East = 1, North, West, South, 
	NorthEast, NorthWest, SouthWest, SouthEast
};


class Cell {
	public:
		double f[Q_DIM];
		double fstream[Q_DIM];
		double feq[Q_DIM];
		bool isObstacle;

		Cell(void);

		double rho(void);
		double rhoOld(void);
		vec<double> u(void);
		vec<double> uOld(void);
		double getFeq(double rho, vec<double> u, int l);
		void setU(vec<double> u, double rho); //poner una velocidad en una celda es poner su poblacion de equilibrio correspondiente

		friend ostream &operator<<(ostream &os, const Cell& cell);
};

class LBM {
	private:
		vec<int> dim;
		double omega;
		double tau;
		double viscosity;
		double soundSpeed;
		vec<bool> periodicBoundary;
	
	public:
		vector< vector<Cell> > grid;
		
		LBM(vec<int> dim, double omega);
		LBM();

		void applyForce();
		void putRectangleObstacle(int x, int y, int dx, int dy);
		void initialize(double rho);
		void stream_step();
		void collide_step();
		double norm();
		double getUMax();

		const vec<int>& getDim() {
			return dim;
		}
		friend ostream &operator<<(ostream &os, const LBM& lbm);
};

#endif
